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i : We compute by numerical transfer-matrix methods the smfaee tree energy r(T). 

Q \ the surface stiffness coefficient k{T), and the single-step free energy s(T) for Ising 

. ^ , ferromagnets with (ooxL) square-lattice and (ooxLxM) cubic-lattice geometries, into 

, which an interface is introduced by imposing antiperiodic or plus/minus boundary 

conditions in one transverse direction. The surface tension a{6,T) per unit length of 
interface projected in the longitudinal direction defines t(T) and k{T) for tempera- 
tures at which the interface is rough: a{9, T)/cose = t{T) + ^K{T)e^ + 0{e^). For 
temperatures at which the interface is smooth, s{T) is the free-energy contribution 
of the dominant entropy-producing fluctuation, a single-step terrace. The finite-size 
scaling behavior of the interfacial correlation length provides the means of investi- 
gating k{T) and s(T). The resulting transfer-matrix estimates are fully consistent 
with previous series and Monte Carlo studies, although current computational tech- 
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nology does not permit transfer-matrix studies of sufficiently large systems to show 
quantitative improvement over the previous estimates. 
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I. INTRODUCTION 



Many phenomena of interest in condensed-matter physics, materials science, and high- 
energy physics occur in systems in which two or more thermodynamic phases coexist below 
some critical temperature Tc. These phenomena are strongly influenced by the behavior of 
interfaces between the phases. Examples of interfacial phenomena from condensed-matter 
physics include wetting, nucleation and growth, and crystal faceting. An example from 
high-energy physics is spatially non-uniform symmetry-breaking in the early universe. A 
number of these systems also have a roughening transition at a temperature Tr (known 
as the roughening temperature), at which the interface changes from microscopically flat 
to microscopically rugged fl],^. For the three-dimensional (3D) cubic-lattice Ising model, 
recent Monte Carlo studies yield T^^ (0.542 ±0.005)Tc ^ 2.45 J/ A;^ The roughening 

transition of one exactly solvable Solid-on- Solid (SOS) model has been shown to belong to 
the Kosterlitz-Thouless (KT) universality class and roughening transitions in related 
SOS models and in the 3D Ising model are generally believed to belong to this universality 
class as well . On the other hand, for the two-dimensional (2D) square- lattice Ising model 
Tr = 0, that is, the interface is always rough. 

An aspect of interfaces that is of great importance in determining interface profiles and 
dynamics is the surface tension, which is the amount by which each unit area of interface 
increases the total free energy of the system in the thermodynamic limit |^. A great deal 
of analytical work , including low-temperature series expansions |[l^ , [T5|JT6[| , has 

been done to study interfaces and their associated surface tensions. The surface tensions 
in a 0"^ field theory [0 and in the g-state Potts model [p!8i|l9| have been studied by Monte 
Carlo simulations. Of particular interest has been the dependence of the surface tension in 
the 2D and 3D Ising models on the orientation of the interface relative to the underlying 
lattice (see Fig. ||); this has been studied both analytically pO|j21|j2^ and through Monte 



Carlo simulations [107,23,24,25,26,27,28 



The analytical treatments of the 3D Ising model rely on approximations which are per- 
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turbative in nature, whereas the non-perturbative Monte Carlo method suffers from critical 
slowing-down and the possibility of becoming trapped in metastable or long-lived unstable 
states. Furthermore, free energies and correlation lengths are cumbersome to calculate in 
Monte Carlo simulations. An alternative method which has been successful in the study 
of phase transitions, particularly in 2D systems, is that of numerical transfer-matrix (TM) 
calculations combined with finite-size-scaling (FSS) analysis [^. The numerical TM FSS 
method does not encounter the aforementioned difficulties: (i) It is nonperturbative and 
does not require foreknowledge of special temperatures (such as Tr and Tc) or of critical 
exponents; (ii) It has no dynamics, so that critical slowing-down does not occur and the 
equilibrium state is readily found; (iii) The free energy corresponds to the largest eigenvalue 
of the TM, and so is readily found; (iv) The calculations are actually easier if the length 
of the system in one direction is infinite; and (v) The correlation length in this direction 
is simple to find since only the two largest eigenvalues of the TM are required. For these 
reasons, there has recently been renewed interest in calculating interfacial properties by TM 
techniques. 

Nevertheless, because the order of the TM increases exponentially with the number of 
sites in a transverse layer (discussed quantitatively in Sec. |I|), the size of systems which 
can be investigated through the numerical diagonalization of a TM is greatly limited. As 
a result, FSS analysis is even more important for TM calculations than for Monte Carlo 
simulations. The order of the TM has also limited most previous attempts at calculating 
the surface-tension anisotropy by TM techniques to capillary- wave and SOS approximations 
2T| , |22|j30| . Specifically, in Refs. pT| , p2| Privman and Svrakic derived FSS relations for the 



correlation length of the Ising model in two and three dimensions by using capillary-wave 
and SOS approximations; these scaling relations involve the surface-tension anisotropy. The 
scaling relations for the 3D model were derived through comparison with the 2D Ising model, 
which has been solved exactly pT| , p2| , p3| , p4| for the boundary conditions used in Refs. pT|j2! 



to generate an interface. In the rough phase, a longitudinal cross-section perpendicular to 
the interface was compared with a 2D system, whereas in the smooth phase the interface 



itself was compared with a 2D system. 

In this article we present a large-scale numerical TM FSS study of surface-tension 
anisotropy in 2D and 3D Ising models. The analytic scaling relations of Privman and 



Svrakic pT] , p2| are applied to analyze free energies and correlation lengths obtained by nu- 
merical diagonalization of transfer matrices describing Ising models on both oox L square 
lattices and oo x L x M cubic lattices into which an interface along the longitudinal 
(x) direction can be introduced by applying appropriate boundary conditions (discussed 
in Sec. 0) in the ?/-direction. The purpose of this study is three-fold: (i) to test the an- 
alytically obtained scaling relations and possibly identify extensions or corrections; (ii) to 
investigate whether the limits on system size imposed by the present generation of parallel 
and vector supercomputers and diagonalization algorithms allow us to reach the asymptotic 
scaling regime in which FSS relations can yield information about interfacial properties in 
the thermodynamic limit; (iii) to compare the resulting estimates of interfacial properties 
with estimates previously obtained by other analytical and numerical methods. Throughout 
this study we use the full Ising model, which includes contributions from microstates with 
bubbles and overhangs as well as SOS microstates. 

We now define the quantities to be studied in this work. Given an interface (see Fig. |ip 
on average perpendicular to the ?/-direction and making only a small angle 6 with the xz- 
plane, the surface tension a{6,T), measured per unit projected area in the x2;-plane, has a 
series expansion for T>Tr [^] 

a{e, r)/cos e = r(T) + -K{T)e^ + o{e^) . (i) 

The coefficients in this expansion are the surface free energy 

r(T) = a(0,T)>0 (2) 

and the surface stiffness coefficient 



k{T) = a{0,T) + 
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(3) 



If T <T{i, K,{T) is infinite, and tlie system can gain entropy only through microscopic fluctu- 
ations, among which single-step terraces are dominant. Each single-step terrace contributes 
a quantity s{T), which is nonzero below Tr, to the total free energy of the system; this 
quantity is known as the step free energy. The step free energy can also be related to the 
surface tension per unit projected area by an expansion when T<Tr [^]: 

ai9,T) /cose = TiT) + siT)\e\ + 0{9^) . (4) 

These three quantities, t(T), k{T), and s(T), are the terms in which the anisotropy of the 
free energy will be discussed. 

The remainder of this paper is organized as follows. We define the model and discuss the 
appropriate boundary conditions in Sec. 0. In Sec. |T| we discuss estimates of the surface 



free energy. Section |^ contains estimates of the surface stiffness coefficient, and Sec. is 
devoted to the step free energy. In Sec. ^ we use the Roomany-Wyld approximant for the 
Callan-Symanzik (3 function p8| , p9| to study the scaling of the correlation length, especially 



near the roughening temperature. Our results are summarized and discussed in Sec. |V11[ 



II. MODEL AND METHODS 

We study the ferromagnetic Ising model with square and cubic lattice geometries. Each 
site i has a spin Sj = ±l which interacts with the spin at each nearest-neighbor site j with 
an interaction constant Ji j > . The Hamiltonian for this system is given by 

Here (i, j) denotes nearest-neighbor pairs on a periodic lattice, and (i) represents boundary 
sites. For nearest-neighbor pairs not on a boundary, Jj = J. Throughout this article the 
lattice constant is of unit length. 

Various boundary conditions are needed to prepare both systems with interfaces and sys- 
tems without interfaces. The presence or absence of an interface is determined by conditions 
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on the boundaries y = l and y = L. All systems we study have the following: (i) a boundary 
interaction constant Jy connecting the site at {x,y = l,z) with the site at {x,y = L,z), (ii) 
a lower boundary magnetic field hi applied on the y = l plane, and (iii) an upper boundary 
magnetic field hi applied on the y = L plane. Table | summarizes the five types of boundary 
conditions used in this paper. 

It should be noted that the "width" of systems with plus/plus and plus/minus boundary 
conditions can be defined by measuring different quantities. These boundary conditions can 
be imposed by coupling the planes at y = 1 and aX y = L to planes of fixed spins. In this 
case there are L+1 couplings in the y-direction, and this number is often used instead of the 



number L of spin sites as the measure of system width (for example, in Refs. [^,^,0]). In 
this article "width" refers to the number of sites L. 

A transfer matrix T can be used to find the free energy per unit length in the longitudinal 
direction. The largest eigenvalue of T, Aq, is related to the partition function Z in the limit 
of infinite longitudinal size (K-^oo) by 0,|T|,|2|,|3|,|9 



hm Z^/^ = lim (Tr 



1/K 



Ao. 



(6) 



In our case the transfer matrix can be chosen to be symmetric, taking into account both the 
interactions within a transverse layer and interactions between adjacent layers. Using the 
thermodynamics of the canonical ensemble, Eq. @ can be used to find quantities such as 
the free energy per site, 

-knT 



f 



A 



In An 



(7) 



where the area A of a layer is given hj A = L in 2D and by A = LM in 3D. 

For the Ising model, a layer in a 2D system consists of a ID lattice with L sites in the 
y-direction, and in 3D a layer is a 2D lattice with L sites in the y-direction and M sites in 
the z-direction. The transfer matrix T is a 2"^ x 2"^ matrix. Exact solutions have been found 
for the free energy and correlation length in finite 2D systems with antiperiodic boundary 



conditions |3^ or plus/minus boundary conditions ||33| , |3^ , so that systems with arbitrarily 



large values of L can be used in FSS. In 3D, on the other hand, numerical diagonalization of 
the large matrix T_ restricts us to small values of A due to limited computer memory. We are 
able to study cross sections as large as A = 19 on a Cray Y-MP/432 by using the Numerical 
Algorithms Group (nag) diagonalization routine F02FJF. The largest cross section we can 
study is y4 = 25, for which an eigenvector of T has 33 554432 elements. Such large systems 



were run on a Thinking Machines Corporation CM2 using a modified power algorithm 
(for details see App. A). 

The longitudinal correlation length is given by by the ratio of the largest and next-largest 
eigenvalues of T_ : 

e = l/ln|Ao/Ai| . (8) 

Other length scales can be defined from smaller eigenvalues of T, and scaling Ansdtze anal- 
ogous to those we use in the next two sections have been proposed [pT| , |22| for these length 



scales. Here we consider only the correlation length defined by Eq. (|^), since it is expected 
that the other length scales require larger values of L (and in 3D, M) to scale in the asymp- 
totic fashion. 



III. THE SURFACE FREE ENERGY 



In order to estimate the surface free energy from finite systems, we compare the free 
energy of two similar systems, one of which has an interface and one of which does not. We 
attribute the difference to the interface. Specifically, we estimate r from antiperiodic (a) 
and periodic (p) boundary conditions in the y-direction by p!0| , pJ]JT^ , |22[ 



Tl{T)^L(fi^\T)-f^>{T)), 



ip), 



(9) 



and from plus/minus (+/— ) and plus/plus (+/+) boundary conditions in the y-direction 

by iigjn|jT|j2| 



rt^-\T)^L{ft^-\T)-fi:/-\T)) 



(10) 



The factor L, which is the width of the system perpendicular to the interface, is needed to 
convert the difference between the volume densities / into a surface density. These relations 
provide just one possible definition of the surface tension in the 2D Ising model, but it 
has been proven |^5| that all such definitions are equivalent with that of Onsager for low 
temperatures. These estimates are valid in any number of dimensions (in 3D the subscript 
L is replaced by L, M). 

Figure ^ shows the 2D estimates t1{T) defined in Eq. (H) and the exact solution t(T) in 
the infinite-lattice limit |^ . The difference between the finite-size estimate and the infinite 
limit in this case is inversely proportional to the correlation length in the periodic system, so 
that for T <Tc, tl(T) converges to t(T) exponentially with L ^2|. Thus the size dependence 



is noticeable only near T^; in all cases, for 2D it has been shown that the critical-point 
behavior is such that <j(T) vanishes linearly with the reduced temperature, t = {Tc — T^jT^. 

The 2D estimates r)^^ \t) defined in Eq. ([T0|) are plotted in Fig. 0b. In this case, 
however, we observe numerically that the estimates from finite-width lattices converge as a 
simple power law: 



(+/-) 



(T) - r(T) ^ A(T)L-^. (11) 



The effective scaling exponent is thus given by 

^) dl^ ' ^^^^ 

where the derivative is calculated by multiple-point finite differences. By holding T = Tc 
and calculating \t) for values of L = 2* x 10, z G {1, 2, . . . , 15}, we can use Aitken's 
method to find lim^^oo |i?(L, T^) — 1| < 10~^, in agreement with the prediction 
limL^oo B{L,Tc) = 1 made in Ref. [^. Using the same method but fixing ksT = 1.5 J, 
we find that for T <Tc, \ limi^^ B{L,T) — 2| < 10^''. For T < and large L, the scaling 
exponent B is observed numerically to be a function only of the critical scaling variable 
tL^/'^=tL (Fig. I). 

The estimates tI^,j{T) for 3D are graphed in Fig. Again the convergence appears 
to be quite rapid except near the Monte Carlo Renormalization Group estimate |^ for Tc. 
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Similar rapid convergence is seen in the estimate Tj^Ii~\t) for 3D (Fig. ^). The value of 
r(T) is exactly known for 3D only at T = Tc (where r = 0) and at T = (where r = 2J), 
but estimates from series expansions improved with Fade approximants [58[ and from Monte 



Carlo simulations [27| agree well with our results for L ^ M. At Tc, r£ appears to converge 
to zero as with little dependence on M, but it is difficult to be certain with only a few 
small systems. Recent analytic results [^] for cubic {K = L = M) systems with antiperiodic 
boundary conditions indicate that the asymptotic scaling region is reached only with very 
large systems. 

IV. THE SURFACE STIFFNESS COEFFICIENT 
A. K in Two Dimensions 



A capillary-wave calculation [0,^] yields that the correlation length for antiperiodic 
boundary conditions and large L should scale as 

for all T < Tc. Solving this equation for k{T) gives a finite-size estimate for the surface 
stiffness coefficient, 

^UT) = ksT-^/. (14) 



Exact results agree with this relation to leading order in L . This estimate is shown 
in Fig. as is the exact solution in the infinite-lattice limit |^2[. We have analyzed the 
convergence and found that for temperatures T < Tc the finite-lattice estimate converges to 
its infinite-lattice limit in a simple power-law fashion: 

kI{T) - k{T) ^ aL-\ (15) 

As in Sec. fTT], we use multiple-point numerical differentiation to determine an effective 
scaling exponent 
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^^^^^^^ jm-niT)] ^ (16) 
d 111 ij 

which is displayed in Fig. |^. Although b is not monotonic, it shares the most important 
characteristics of the scaling exponent B defined in Eq. ([T2|); that is, b varies from 6=1 for 
T = Tc (in agreement with critical scaling) to 6 = 2 for T'^T^, and 6 is a function only of the 
critical scaling variable tL for large L. 

The exact solution for [^2(T)]~^ is a sum in which both the total number of terms 



and each individual term depends on L [^. We avoid the resulting difficulty of exactly 



differentiating with respect to L by using numerical differentiation in the evaluation of 



Eq. ([Tq). However, it is also possible to define a "reduced correlation length" ^, obtained 



from the antiperiodic correlation length and the periodic correlation length by 



which contains only one explicitly L-dependent term (see Eq. ( pTf ) in App. B). Replacing 
^£(T) in Eq. (|T^) with ^l(T), we find an estimate, KLiT), which is a differentiable function of 
ln(L) if L is treated as a continuous variable. Except for small systems or T^Tc, [C,l{T)]~^ 
is vanishingly small in comparison with [^£(T)]^^. Thus it is clear that ^^{T) will be less 
than kIC^) but must converge to the same limit. In fact, for all L, k,{T) < kl(T) < k,1{T). 
In the limit of large L we find the effective scaling exponent 

/ 1 \ 



b{tL) 



4c(tL)2[^(tL) -2c]/7r2 

2 



TT 



2{tL)-9{tLMtL)-2c]) ' ^^^^^ 



where c = sinh ^(1) ^ 0.8814 and g{x) is given by 



7r^2 



X 



g{x) = d{2cy+ -) . (18b) 



Equation ( |18a| ) is derived in App. B. As is the case for B(tL) and b{tL), b{tL) varies from 



1 at tL = to 2 at tL = oo, but unlike b{tL), b{tL) is monotonic (Fig. ^d). It is proved 
in App. B that \imtL^oob(tL) = 2, and since limtL_+oo 6(tL) = \imtL-^oob{tL), we also have 
lim^i^oo b(tL) = 2. 
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Another capillary-wave calculation pl|j2^ indicates that for plus/minus boundary con- 
ditions the correlation length scales as 

for all T <Tc, giving the finite-size estimate for the surface stiffness coefficient 

4"/-^(r) = fc.T^^^. (20) 

Figure shows k,^^^ We have observed numerically that for any fixed temperature 

T<Tc, k^]^^^\t) converges to k(T) as L~^. 

Refs. and p2| actually state Eqs. (|19D for large L with ^ ''^ (x . However, in the 
same papers Privman and Svrakic obtained estimates for the correlation length in 3D both 
for T > Tr (Eq. (p2D) and for T < Tr (Eq. (p6D), which should be related to the scaling 
relation of the correlation length in 2D (Subsec. [IV B|) . It is reasonable to suspect that the 
L-dependence of the correlation length is not affected by the roughening transition in the 
case of plus/minus boundary conditions, just as it is not affected in the case of antiperiodic 
boundary conditions (Eqs. (|2T| ) and (^)). This argument and our numerical observations 
lead us to conclude ^i^^""* oc (L + 2)^. 



B. n in Three Dimensions 



There are two contributions to the correlation length along the interface. The first comes 
from pieces of interface which, due to thermal fluctuations, are normal to the longitudinal 
direction, which is the direction in which the correlation length is measured. This con- 
tribution is related to the surface stiffness k(T). The second contribution comes from the 
formation of bubbles in the bulk and is much less important for T^T^. 

Both of these contributions can be treated approximately by considering a cross section 
of the 3D system parallel to the xy-plane. The relationship between correlation length and 
surface stiffness in such a cross-section should be the same as for a 2D system of the same 
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width L. Since there are M such cross-sections, and since each contributes on average the 
same amount of interface normal to the longitudinal direction, the total effective surface 
stiffness for the 3D system is given by Mk{T) and has the same form as the corresponding 
2D systems; that is, the appropriate estimate for the 3D surface stiffness coefficient for 
antiperiodic boundary conditions is given from Eq. (|T^ by 

^Im{T) = ^BTy , (21) 



and for plus/minus boundary conditions from Eq. ( pOD by 

^L,M (T) - ksT—jj-—^ , (22) 

for Lk. M^l. These estimates were proposed by Privman and Svrakic |21|j2^. No part of 
this argument depends on the boundary conditions in the z-direction, so Eqs. (pi]) and (^2]) 



are valid for either periodic or free boundary conditions applied in the z-direction. In this 
article we discuss k(T) only in the case of periodic boundary conditions in the z-direction, 
but we have observed that little quantitative change occurs when free boundary conditions 
are used. 

In Fig. ^ the "reduced correlation length," defined in 3D by analogy with Eq. (0), is 
used in place of ^£ mO^) estimate k(T) using Eq. (|2ll) . As in Subsec. [IV7\| , this produces 



a reduction of the estimate of n near (where k(T) is zero in an infinite system) but cannot 
have any effect on the estimate of k(T) for T <Tc in the limit of large systems. Also shown 
in Fig. 1^ are estimates of k(T) based on series expansions improved by Fade approximants 
48| and estimates based on Monte Carlo simulations In spite of the small systems to 



which our study was limited, we find good agreement between the transfer-matrix estimates 



in this study and the previous estimates |^8|,^ over a wide range of temperatures between 
Tr and Tc. 

In Fig. 1^ we show the estimate for k{T) that results from using Eq. (|22D. Comparison 
with the series expansion and Monte Carlo estimates shows that finite-size corrections to 
Eq. ( p2|) cannot be neglected at any temperature for the systems we studied. 

13 



The surface stiffness coefficient for the infinite system shows a characteristic KT disconti- 
nuity at Tr, jumping from k{T) = oo for all T <Tr to K{T)/kBT = 7t/2 for T = T^ [jST] . 
If we use the criterion h{Tr) /kBTu = 7r/2 to estimate Tr while using antiperiodic boundary 
conditions, we find that the estimated value of Tr depends most strongly on the value of 
M (Fig. For the system {L,M) = (3,8), we find Tr ~ 2.337 J/ks- This represents a 
deviation of approximately 5% from the Monte Carlo estimates [0,^,^, and is comparable 



with our alternative estimates of Tr discussed in Sec. Although the convergence in Fig. § 
appears to be roughly power-law in behavior, much slower convergence should be expected 



for large systems due to the KT nature of the roughening transition [52 



V. THE STEP FREE ENERGY 



For T < Tr, a SOS calculation indicates that the correlation length for antiperiodic 



boundary conditions should scale as pl|j2^ 



eLM{T)^fi{T)M-exp 



Ms{T) 



TT ^ -2 



(sm-l , (23) 



ksT 

where /i(T) = ^ exp[wC(T, L/M)], and C{T, L/M) is a function possibly of the temperature, 
boundary conditions, and shape of the layer, but is independent of the size of the layer. The 
boundary conditions in the 2;-direction may be either periodic (Fig. |^a) or free (Fig. ^b). By 
taking three different values of M we can solve simultaneously for s(T), wlT), and fi{T). 
The step free energy s(T) discussed in this article is the step free energy for an interface 
parallel to the x2;-plane. A low-temperature expansion for the step free energy of inclined 
planes can be found in Ref. [T5|| . 

Unfortunately, numerical convergence difficulties limit to very small systems our studies 
of interfaces at very low temperatures. Furthermore, quantities near Tr, where we can use 
larger systems, may be expected to be difficult to estimate due to the KT transition. Nev- 
ertheless, the estimates of s{T) rapidly converge at low temperatures, so that the limitation 
to small systems may not present difficulties in estimating the step free energy. Also, the 
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estimates of s{T) are rather insensitive to the values of w and fi{T), so that problems at low 
temperatures with the last two quantities (described below) do not seem to create difficul- 
ties in the estimates of s{T). This is because the divergence of the correlation length at low 
temperatures is dominated by the exponential factor in Eq. p3|), which contains the step 
free energy. 

The estimates of s{T) derived from Eq. (|23|) are compared in Fig. ^with approximations 
valid either near T = or near T = Tr. At low temperatures the interface becomes a series 
of plateaus separated in height from their neighbors by one lattice spacing. These plateaus 



are equivalent to domains in a 2D system |15,22|; the equivalence is exact at T = 0. This 



yields the well-known approximation |15 



s(T) ^ r2D(T) . (24) 

We use r^^(T) as an approximation to s(T) for T <^Tji for the following reasons: it provides 
a simple physical picture; it agrees with the first several terms of low-temperature series 
expansions of s{T) ||T5|; it has been solved exactly |31|; and the 2D critical temperature 



is a rigorous lower bound on Tr [^. At higher temperatures (T ^ Tr), s{T) shows the 
characteristic KT essential singularity at Tr PJ^H], and is given by 



where the nonuniversal parameter c has been found by Monte Carlo simulations [0,3 to 
be c=1.57± 0.07. In Fig. ^we have calculated the estimate Eq. ( ^5]) with c = tx/2. 

The correspondence between a low-temperature interface in 3D and a 2D system ensures 
|55| , |56|j57| , |58[| that for periodic boundary conditions in the 2;-direction w = 1/2, whereas for 
free boundary conditions in the 2;-direction w = Q. Our estimates, however, have w = Q for 
sufficiently low temperatures for both sets of boundary conditions (Fig. |I0|). The incorrect 
estimate for w in the case of periodic boundary conditions in the z-direction provides a 
caution that the asymptotic scaling regime, in which corrections to Eq. ( p3D can be expected 
to be very small, has not been reached with the current small systems. The parameter /i(T) 
(not shown) is consistent with w = at low temperatures. 
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For plus/minus boundary conditions, the scaling relation 



eW/"^(T)^^(r)M-exp 



Ms(T) 



sm — r sm 



IT 



1 -1 



2(L + 2) 2(L + 2) 



(26) 



is expected in analogy with Eq. (^) on the basis of another SOS calculation |2l|j2^ . The 



analysis for Eq. ( P6D parallels that for Eq. (|23|), so that at zero temperature we should find 
that for periodic boundary conditions in the 2r-direction w = l/2, and that for free boundary 
conditions w = 0. In this case the values extracted for s{T) (not shown) for different system 
sizes are inconsistent with each other even at low temperatures, and the corresponding values 
of w (also not shown) show no pattern whatsoever. 

In short, the scaling relation of Eq. (pSf ) yields values of the step free energy that are 
consistent with the low-temperature approximation r^°(T), whereas Eq. (^4]) does not even 
yield consistent estimates at low temperatures. The parameter w{T) is found to be wrong 
at low temperatures when calculated from Eq. ( ^3] ) with periodic boundary conditions in 
the 2;-direction, and when calculated from Eq. (|2^ ) it yields inconsistent values. It is thus 
apparent that large corrections to scaling are present for the systems in this study. 



VI. THE CALLAN-SYMANZIK /3 FUNCTION IN THE ROUGH PHASE 



In Subsec. [IV B| and Sec. we have evaluated terms in the surface tension from FSS 
relations for the correlation length |^I| , p2| in 3D for 0< TkT^. As discussed in Subsec. [IV B 



for Tji<T <Tc each of the M cross-sections parallel with the xy-plane makes a contribution 
to the correlation length in an ooxLxM system which is of the same form as the correlation 
length in a 2D system of width L. From this argument we have Eqs. ([2T[) and (p2[), which 
give ^L,M oc M for either antiperiodic or plus/minus boundary conditions. Such linear 



divergence with system size of the correlation length is characteristic of a critical point [29 



and indicates that the rough phase of the interface is a critical phase [0. A quantity that 
has proven particularly useful in locating such extended critical phases in previous studies 
{e.g. Refs. |BP[,KD[) of systems with Kosterlitz-Thouless behavior is the Roomany-Wyld 



approximant to the Callan-Symanzik (3 function [^,^. This approximant is defined by 
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^"'^^^^ " {dHi,,M/L-])/{dT) ' ^^^^ 
where the derivative with respect to InM is evaluated by multiple-point finite differences. 

(There should be no confusion between I3l,m{T) and the dimensionless inverse temperature, 
which is called (3 in the appendices.) Within the region where a "critical interface" exists, 
one should have I3l,m{T) = const ^ 0. In this section we apply the jS function to the 3D 
Ising model with antiperiodic boundary conditions in order to obtain independent TM FSS 
estimates both for Tr and for the (size- dependent) temperature where the rough-interface 
phase is destroyed by critical-point fiuctuations in the small systems we are studying. 

In Fig. |ll| we show various estimates of I3l,m{T) for the case of antiperiodic boundary 
conditions in the y-direction and periodic boundary conditions in the z-direction. Finite- 
size effects were larger for the case of free boundary conditions in the ^-direction and much 
larger for the case of plus/minus boundary conditions in the ?/-direction, so these results are 
not shown. The dashed curve near T = is the value of lim^^^M^oo /^L.Af (7") where r^^(T) is 
substituted for s{T) in Eq. (^3]). If the scaling relations of Eqs. (^) and (|26| ) for T <Tr 
are to join continuously with the scaling relations of Eqs. (^) and (|22| ) for T >Tr at Tr, 
then w(Tr) = 1. The solid curve for T<Tr traces ImiL^M^oo I3l,mO^) under the assumptions 
that Eq. (^31) holds near Tr with w = 1, and that the step free energy near Tr is given by 
Eq. (p5|) . The agreement of the data in Fig. |ll] with this portion of the solid curve is good 



despite the fact that Eq. (p3|) is based on the assumption that the terraces in the interface 
are well separated, which is not the case near Tr. The solid horizontal line I3l,m{T) = for 
the temperature range Tr<T <Tc shows liniL^M^oo /3l,m(7') assuming that Eq. (^If ) holds 
in this range. 

The finite-size estimates for PlmC^) °f types. The first type consists of evalu- 
ating ((iln[^2,_Af/L^])/((ilnM) by a three-point estimate while holding L constant, and then 
making the empirical extrapolation 

/3l,(m^oo){T) = (3lm{T) + a(T) exp[-6(r)M] (28) 

where a(T) and 6(T) are unknowns that depend only on temperature. (The only justifi- 
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cation for this extrapolation procedure is that it appears to work nicely for small systems, 
although we expect logarithmic convergence for large systems because roughening is a KT 
transition [Q.) For L>3 this yields results which depend little on M and show the emer- 
gence of a plateau near Tr with i3l,m{T) ^ —0.06. Even for L = 5, critical finite-L effects 
destroy the plateau for temperatures greater than about 3.2 J/ ks- The extrapolated value 
of Pl,{m^oo)(T) for L = 3 does not show any plateau, but there is a change in concavity at 
about 2AlJ/kB, where the critical finite-L effect gives way to the T kTr behavior. 

The second type of finite-size estimate is made by using only systems with square cross- 
sections (L = M) in the numerical evaluation of ((iln[^ijv//L^])/((ilnM) in Eq. (|27|) . This 
estimate of Pl,m(T) yields a plateau at Pl^m{T)^—0.1. 

In either case, the onset of the plateau (which corresponds to an estimate of Tr) can be 
defined by a peak in dP [3L^MiT) / dT"^ . This estimate yields Tr ^ 2.3 J/ ks, which is approx- 



imately the same as that given in Subsec. [IV B| . The high-temperature end of the plateau 
seems to correspond to the temperature at which critical finite-L effects become large in the 
interface internal energy (not shown), which is given by 



.2diT/T) 



U = -T'^^ . (29) 



VII. DISCUSSION AND CONCLUSION 

As a nonperturbative numerical method for studying semi-infinite systems in equilibrium, 
numerical transfer-matrix (TM) calculations followed by finite-size scaling (FSS) pl|j2^ 
provide an appealing alternative to the analytical and Monte Carlo methods which have 
been used to study interfaces in the Ising model. Nevertheless, the rapid growth of the order 
of the transfer matrix with system size has limited most previous studies to Solid-on-Solid 
(SOS) and capillary-wave approximations. 

In this paper we have performed large-scale numerical TM calculations on a Cray Y- 
MP/432 and a Thinking Machines Corporation CM2 to obtain the correlation lengths and 
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free energies of 2D and 3D Ising ferromagnets with interfaces, which are related to the surface 



tension anisotropy through the scahng relations derived in Refs [^,^]. In 2D we are able 
to analyze the corrections to scaling. In 3D we are only able to confirm the scaling relations 
for antiperiodic boundary conditions, and the resulting estimates for t{T) and /t(T) agree 



well with estimates from series expansions improved with Fade approximants and from 
Monte Carlo simulations |2^. For plus/minus boundary conditions the finite-size effects in 



our estimates are still large, even for the largest systems which can be studied using current 
supercomputers and algorithms, and the asymptotic scaling regime has not been reached. 

We find that although in 2D the surface free energy converges exponentially for an- 
tiperiodic boundary conditions, for plus/minus boundary conditions the convergence obeys 
a power law. The effective scaling exponent for plus/minus boundary conditions is a func- 
tion only of the critical scaling variable tL, and varies monotonically from 1 at tL = to 2 as 
tL^oo. Although the exact value of t(T) is not known in 3D for arbitrary temperatures, 
our numerical results for L = M agree well with the estimates of Refs. p^j48[1 . We are not 
able to make a quantitative study of the scaling. 

In 2D the surface stiffness coefficient is found from scaling relations based on capillary- 
wave studies pT| , |2^ . For antiperiodic boundary conditions we find that n1{T) at any fixed 



temperature T < Tc converges in a power-law fashion with an effective exponent similar 
to that found for r^+/^)(T). Specifically, the exponent is a function of the critical scaling 
variable tL, and varies from 1 at tL = to 2 as — oo, though not monotonically. By 
using the "reduced" correlation length, a monotonic exponent, discussed in App. B, can 
be calculated in the infinite-lattice limit. For plus/minus boundary conditions Ki^^~\t) 
converges simply as 1/L. 

The estimates for n{T) in 3D are based on scaling relations similar to those in 2D because 



of the similarity of a cross-section in 3D to a 2D system. These estimates pT],p2| in the case 



of antiperiodic boundary conditions agree well with both series and Monte Carlo ||27 
estimates. Plus/minus boundary conditions lead to large finite-size effects in the estimates 
of k(T) at all temperatures. Estimates for the roughening temperature based on the KT 
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relation n{Tii) /ksT^ = 7r/2 |]7|j8|j5T[] depend more strongly on the value of M than on the 
value of L. For (L, M) = (3,8), this estimate yields Tr^2.3. 



Solid-on-Solid arguments have yielded scaling relations for the correlation length pl|j22 
which we exploit to find s{T). In the case of antiperiodic boundary conditions, the estimates 
for s(T) agree well with other approximations at low temperatures. The low-temperature 
value of the fitting parameter w for systems with periodic boundary conditions in the z- 
direction, however, is in contradiction with known results |l55| , |56| , |57| , [58[| . This indicates the 
presence of corrections to scaling for the correlation length in the small systems of this study. 
These corrections are even more pronounced in the case of plus/minus boundary conditions, 
where the values of s{T) estimated from different system sizes are inconsistent even at low 
temperatures. 

The rough phase of the interface, which is a type of critical phase [ffl], displays scaling 



behavior for the correlation length which is similar to that of other critical systems [29 



Specifically, as a result of Eq. ( |21| ) , we have ^ oc M. We used the Roomany- Wyld approximant 
to search for a plateau in the Callan-Symanzik (3 function PB| , |55| which corresponds to this 
critical temperature range. The onset of the plateau agrees roughly with the estimate of 
Tr from Subsec. p.V lj| , while the end of the plateau corresponds to the onset of large-scale 
finite-size corrections in the interface internal energy caused by the disappearance of the 
interface at T^. 

In summary, we have provided transfer-matrix finite-size scaling estimates of the quan- 
tities which characterize the surface-tension anisotropy in two- and three-dimensional Ising 
ferromagnets: the surface free energy r(T), the surface stiffness coefficient k(T), and the 



step free energy s(T). Our results support the analytically obtained relations of Refs. pl|j22 



and are fully consistent with existing series and Monte Carlo estimates [^,^]. However, we 
find that in the three-dimensional case the limitations of current computer technology and 
algorithms do not permit transfer-matrix studies of sufficiently large systems to improve on 
the numerical accuracy of the previously existing estimates. 



20 



ACKNOWLEDGMENTS 



The authors wish to thank K. K. Mon and V. Privman for useful discussions, M. Grant 
for a reading of the manuscript, and P. E. Oppenheimer for discussions about PARIS pro- 
gramming of the CM-2. The authors would also like to thank M. E. Fisher for assistance 
in locating relevant research and for providing unpublished independent series estimates 
from his work with H. Wen, and would like to thank M. Hasenbusch and K. Pinn for pro- 
viding unpublished independent Monte Carlo estimates for comparison. This research was 
supported in part by the Florida State University Supercomputer Computations Research 
Institute which is partially funded by the U.S. Department of Energy through contract ^ 
DE-FC05-85ER25000, by the FSU Center for Materials Research and Technology, and by 
the National Science Foundation through grant # DMR-9013107. 



This appendix contains details about the algorithms used in the diagonalization and the 
programming implementation of these algorithms. In this appendix we take kB = J= 1, so 
that all quantities are dimensionless. 



can be decomposed into the product of two matrices, T = AD . The matrix D., which can 
be chosen to be diagonal, takes into account both interactions between spins in the same 
layer and interactions of the spins with a magnetic field. The matrix A takes into account 
the interaction between one layer and the next layer. It is a direct (Kronecker) product of 
A identical 2x2 matrices 



where j3 = l/T. In both of the algorithms used, as described below, the core of the algorithm 
demands that for any given vector v one can obtain a vector w = Tv. In order to be able 
to find even a few of the largest eigenvalues of T with A large, one needs to minimize the 



APPENDIX A: COMPUTATIONAL DETAILS 



If each of the identical layers has A spin sites, the 2"^ x 2"^ transfer matrix T in Eq. 




(Al) 
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computer storage required. This means that only a few (typically between 2 and 20) vectors 
of size 1^ can be stored, whereas the entire 2^^ elements of T cannot be stored. The matrix 
D_ only requires storage of 1^ elements, so it is computed once for each temperature and 
stored. All of the boundary conditions used in this paper are incorporated into the matrix 
D_. Two slightly different methods of computing Ax were used. 

Although the size of the matrix T can be reduced by consideration of the symmetries of 
the model, the savings in computer storage are not great in systems into which an interface, 
which breaks many symmetries, has been introduced. Furthermore, the symmetry-reduced 
matrix would reduce the speed of the diagonalization program discussed below on the CM- 
2 by complicating interprocessor communications. Since very few systems with reasonable 
aspect ratios could be added by using symmetry reduction, and because it would be inefficient 
to implement symmetry reduction on the systems now accessible, we have not pursued this. 

For A< 19, the NAG routine F02FJF was used on the FSU YMP, which has 32 MWords 
of memory. This routine required use of the symmetric transfer matrix D^^^AD^^^- One 
can write A = Aj^Aj^_i ■ ■ ■ Ai, where each 2^ x 2"^ matrix A^ is a direct product of A — 1 
identity matrices of size 2x2 and the matrix a in position i. A permutation matrix P_ can 
be defined such that P.^^Ai P_ = A^^i with Ay^^i = A^. In other words, P permutes the 



matrices within the Kronecker product in a cyclic fashion [|61i|64|6^1 |64| , |65 i |66| . For any j, it 
follows that 4. = P'^^^A^P^-\ and thus one obtains A = p-'^^^A^ (PA^)^-\ Since P'^ 
is the identity matrix, all one needs to do is to program the multiplication of a vector by 
Ai and by P and repeat these A times to obtain the multiplication by A. The matrix Ai 
has only two non-zero elements in each row and column, so multiplication by A-^ can be 
efficiently programmed. One way of doing this is to break the vector to be multiplied into 
an upper and a lower part, each having 2^~^ elements. Then 

\VbJ \Vb/ \Va) 

The matrix P acts on a vector as = with the iih. even element of w given by 

the ith element of and the ith odd element of w given by the iih element of Va- The 
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multiplication by P_ can be efficiently programmed using scatter-gather routines on a vector 
computer. 

For 14 < y4 < 25, a double-precision program in FORTRAN/PARIS for a Thinking 
Machines Corporation massively parallel CM-2 computer was used for the diagonalization. 
The use of the PARIS programming language was essential for diagonalization of the large 
matrices, since the FORTRAN compiler would have otherwise allocated more than the 
minimum number of arrays. Programming simplicity in the PARIS portion of the program 
required that the number of 1-bit processors used be less than or equal to the size of a single 
vector. The CM-2 available has 2^^ processors, and can be divided into quarters with 2^^ 
processors each, giving the limit 14 < A. The upper limit was imposed by the computer 
memory on the CM-2. To calculate the correlation length, at least two vectors need to 
be iterated. This entailed storing 6 double-precision (64-bit) vectors: one of these a work 
vector, one the vector associated with the diagonal matrix D_, and two vectors each for 
the current iteration vector and for the next iteration vector. In addition, A — 1 single-bit 
vectors were used to efficiently implement the multiplication of a vector by the matrices A^. 
(Temporary vectors used in the formation of the bit vectors and of the diagonal matrix were 
discarded before the 64-bit iteration vectors were allocated.) For A = 25, this required use of 
1.7 Gbytes of the available 2 Gbytes of main memory on the CM-2. The implementation of 
the multiplication by one of the A matrices A^ was made using the hypercube communication 
of the CM-2. This was done by using the decomposition A^ = e^J -|- e~^Xj. Here / is the 
2'^x2'^ identity matrix, and the permutation matrix 2Li is a direct product of A—1 identity 
matrices of size 2x2 and the matrix X = (J^ in the ith position. The multiplication by 
the 2Li matrices was performed using efficient communication with the PARIS command 
GET_FROM_POWER_TWO, which allows hypercube communication to elements a power 
of two away. 

The program was written to diagonalize any real matrix with both the largest eigenvalues 
and corresponding eigenvectors real. The algorithm used is based on a generalized power 
method and proceeds as follows [^. A 2^xR matrix U_ is initialized with random elements. 
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Due to constraints of the PARIS program, R must be an integer power of 2, and memory 
constraints dictate that R = 2 if A = 25. Then the matrix V_ = TU is obtained using 

the muhiphcation procedure described above for each of the R vectors in V_. Two Rx R 
matrices S = UjH and Q = U^V are then formed. The eigenvalue equation S_ ^ QG = QT> 
is then solved to find the RxR diagonal matrix of eigenvalues V and the RxR orthogonal 
matrix of normalized eigenvectors Q_. The diagonal elements of V provide the estimates for 
the eigenvalues of the transfer matrix T for this iteration. The matrix U_ is then updated by 
U_ ^ VQ , and the process is repeated. 

Near Tc typically about 25 iterations were necessary to obtain convergence of the corre- 
lation length to one part in 10^^. However, at lower temperatures many more iterations were 
required for convergence. Near Tc for A = 25, the computer time required for convergence on 
the 2^^ processor CM-2 with double-precision Weitek fioating-point accelerators was about 
15 minutes. 



APPENDIX B: THE SCALING EXPONENT FOR k IN 2D 



In Eq. (|18a|) in Subsec. (|IV A|) we presented the exact result for an effective scaling 
exponent b{L, T) such that 

kL{T) - k{T) ^ aL-^ . (Bl) 

Then, in analogy with Eq. ([16|) , 

6(L,T)^^^MMnZ^ (B2) 
dmL 



kL{T) - k{T) dL 



This appendix contains a sketch of a derivation of Eq. ( |18a| ), which is made by rewriting 
the functions in Eq. ( [B3| ) in terms of the reduced temperature t = (Tc — T) /T^ and the 
critical scaling variable x = tL, and then keeping only the 0(1) term in t at fixed x. In this 
appendix we take kB = J=l and measure the surface stiffness coefficient in units of ksT, so 
that all quantities are dimensionless. 
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First, we define the functions wliicli are to be expanded. Tfiese are pT| , |22|j3^ tlie inverse 
temperature 



/? = 1/T , (B4) 



tlie Onsager angles 



7(0) = 2(3 + ln(tanli(/5)) (B5a) 



and 



7(7r/L) = coshr\A - cos(7r/L)) , (B5b) 

wliere 

^ cosIi^(2/j) 

^ - sinli(2/3) ' ^^^^ 

tlie "reduced" mass gap, which is the sum of the mass gaps for antiperiodic and periodic 
boundary conditions, 

= 7(7r/L)-7(0), (B7) 
the reduced surface stiffness coefficient 

and the surface stiffness coefficient in the limit L — 00 

K = sinh7(0). (B9) 



Note that although L was initially allowed to take only integer values in Eqs. ( [B5b| )-( P8D , 
we drop this constraint to find a function /«l(T) which is smooth and differentiable with 
respect to L. 

Inserting Eqs. (|B|-^) into Eq. (|B3D, we find 
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dkL{T) _ r-2 1 rfmz, 1 

dL 2L?rhL \ L tjil dL ] 
. f-2 1 d7(7r/L)l 
L rhi dL J 



-2k 



-2 TT ci7(7r/L) 1 
T ^ L2mL d{7i/L) J 

sin(7r/L) 



1 
L 



^ - cos(7r/L))2 - 1 



so that 



b{L,T) 



2LkL{T) 1 kL 



sin(7r/L) 



/€l(T) -k(T) vr ^(^ _ cos(7r/L))2 - 1 

It is convenient to define 

2 



sinh ^(1 



A little algebra then yields the Taylor expansion 

A ^2 + 2cH'^ . 

With this result, and with L = x/t, it is easy to find that 

cosh7(7r/L) = A — cos ( - 

^ (2 + 2ch'^) - ( 1 



■nt 

X 



n 2^ 



1 + -g\x)t\ 



where for convenience we have defined 

g{x) = \/(2c)2+Q' 
Now taking cosh^^ of Eq. (|B13| ) gives 



.-1 



(^A — cos 








) 




. X . 





g{x)t . 
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We also expand 7(0) in a Taylor series to give 



7(0) = 2/3 + In (tanh[/3]) 



7(0) 
+ 2ct . 



^7(0) 
t=o «t 



t=o 



From Eq. (|B9| ) this implies 



Using Eqs. (|B7D, O, (pT5| ), and pT6| ) provides 



2x2{7(vr/L) -7(0)} 



2x'^{g{x) - 2c} 
In order to evaluate b{x) it is also necessary to utilize 



dl{n/L) 
di-K/L) 



sin(7r/L) 



[A - cos{n /L)Y - 1 
(vrt/x) 



7C 



xg{x) ■ 



Then, using Eq. ( B10a| ) 
dL 



-2 71 dj{n/L)' 
T ^ L'^rh d{7i/L) 



-2t 2 . TT 



2x'^[g{x) — 2c] [ X Tc xg{x) 
Using this and Eq. (|B18| ) we find 



dL 



-2t 2kL 
+ 



2x'^[g{x) — 2c\ \ X xg{x) 



2x^[g{x) - 2c] 



X 1 



TC 



2x'^g{x)[g{x) — 2c] 
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Finally we can combine Eqs. (0), (|BT7| ), ( pTHD , and ( p2TD to obtain 

L dkLiT) 



b{x) 



kL{T) - k(T) dL 
( 7iH/2x^[g{x) -2c] 



\{7TH/2x^[g{x) - 2c]} - 2ct 

V 2a;2(7(x)[(^(a;)-2c]; 

^ (l -4cx2[^(x) -2c]/7r2) 
X 1- 



TT 



(B22) 



2a;^(yf(a;)[(yf(a;) — 2c] y ' 

which is correct to 0(1) in t. Since we are interested in b{x) in the limit t — > (which for 
fixed X necessarily means L oo), this is the final answer. 

It is instructive to determine b{x) for x = and x oo. Since we can write 



xg{x) = ^ {2cxY + 77^ 



(B23) 



we immediately find 



lim xgix) = it. 
> 



Substituting this into Eq. (|B22|), we find 



(B24) 



b{x = 0) = 1, 



(B25) 



which agrees with critical scaling p9[. We also note that 



lim x'^\q(x) — 2c] = lim c ( — ] 



1 - 



1 / 71 



4 \2cx 



so that Eq. ( p22D yields 



(B26) 



lim b{x) = 2. 



(B27) 
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FIGURES 

FIG. 1. A 3D system illustrating the axes and dimensions of the system, and the orientation 
of the interface, which separates the region of positive magnetization from the region of negative 
magnetization. In this figure boundary conditions are used to tilt the interface by a fixed angle 9 
in the z-direction. Although the boundary conditions discussed in the text do not produce such a 
tilt, they can be used to predict the change a tilt would produce in the surface tension. 

FIG. 2. The estimates of the surface free energy in two dimensions, (a) r^(r), from Eq. @, 
and (b) t^j^''\t), from Eq. @. The solid curve is the limit of ti (T) for an infinite lattice ||3l| . 
The solid vertical line gives the 2D critical temperature. 

FIG. 3. The effective scaling exponent B(tL) for Tj^^ \t) in two dimensions, defined in 
Eq. (|l2|) and evaluated by numerical differentiation. 



FIG. 4. The estimates of the surface free energy in three dimensions, (a) t£ from Eq. 

and (b) t^^j from Eq. (p!^). The dashed curve is a low-temperature expansion the solid 

curve comes from series expansion improved with Fade approximants E^, and small solid circles 



(bigger than the error bars) indicate the results of Monte Carlo simulations |27]. The solid vertical 
line marks the 3D critical temperature as determined by a Monte Carlo Renormalization Group 
study @. 

FIG. 5. The estimates of the surface stiffness coefficient in two dimensions, (a) K,'j^{T), from 
Eq. (0), and (b) from Eq. (pO[). The solid curve is the limit of kl(T) for an infinite 

lattice. The solid vertical line is the 2D critical temperature. 

FIG. 6. The effective scaling exponents for the surface stiffness coefficient in two dimensions. 



(a) The scaling exponent b{tL), defined in Eq. (16), for k'^{T). (b) The scaling exponent b{tL), 



defined in Eq. (|B3|), for kL{T). The solid curve is the limit for large L referenced in Eq. ( |18aD and 
derived in App. B. 
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FIG. 7. The estimates of the surface stiffness coefficient in three dimensions, (a) kL,M{T), from 
Eq. (21), and (b) from Eq. (^). The solid curve comes from series expansions improved 



with Fade approximants |4^. The small solid circles (bigger than the error bars) are the results of 
Monte Carlo simulations |27]. The vertical dotted line indicates the roughening temperature, Tr, 



as determined from Monte Carlo studies [Q]. The vertical solid line is the 3D critical temperature 

FIG. 8. The fractional difference between the transfer- matrix roughening temperature, T^^, 
as determined from KL,MiT)/kBTji = tt/2, and the roughening temperature found from by Monte 
Carlo simulations Q, T^^. The error bars (shown only for L = 3) result from uncertainty in T^'~'. 

FIG. 9. The estimates for the step free energy using antiperiodic boundary conditions in the 
y-direction. Estimates using periodic and free boundary conditions in the z-direction are given by 
(a) s'*P(T) and (b) s^^{T), respectively. By keeping the temperature constant and taking 3 different 
values of M (and any combination of values of L), we solve Eq. ( |23D simultaneously for s, w, and 
fj.. The solid curve is the 2D surface free energy. The dashed curve shows the estimate Eq. (25) 
|23,|4|. 



FIG. 10. The estimates for the exponent w, found by solving Eq. (|2!^ ) for s(T), w, and ^(T), 
where the boundary conditions in the y-direction are antiperiodic. It is known that at T = 
■u; = l/2 for periodic boundary conditions (ap) in the z-direction and that 'W = for free boundary 
conditions (af). At low temperatures our transfer-matrix estimates give that w = for both sets of 
boundary conditions, indicating the presence of substantial corrections to scaling. 



FIG. 11. Finite-size estimates for the Callan-Symanzik /3 function, given in Eq. (|27|). Results 
are shown both for sets of systems with square cross sections and for extrapolations based on 



Eq. (|2q). The solid curve comes from inserting the estimate of Eq. (25) |£4| for s{T) for T^Tji 



into Eq. (p3). The dashed curve comes from an approximation of the step free energy by the 2D 



surface free energy in Eq. (23). 
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TABLES 



TABLE I. Boundary conditions discussed in the text. 



Name Jy hi Hl 

Free 

Periodic J 

Antiperiodic —J 

Plus/Plus J J 

Plus/Minus J -J 
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